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p^ Abstract 

O, We study a trajectory-planning problem whose solution path evolves by means of 

<^ a Lie group action and passes near a designated set of target positions at particular 

^ times. This is a higher-order variational problem in optimal control, motivated by 

,—1 potential applications in computational anatomy and quantum control. Reduction by 

symmetry in such problems naturally summons methods from Lie group theory and 

l/^ Riemannian geometry. A geometrically illuminating form of the Euler-Lagrange equa- 

Q tions is obtained from a higher-order Hamilton-Pontryagin variational formulation. In 

this context, the previously known node equations are recovered with a new interpre- 
tation as Legendre-Ostrogradsky momenta possessing certain conservation properties. 
Three example applications are discussed as well as a numerical integration scheme that 
I I follows naturally from the Hamilton-Pontryagin principle and preserves the geometric 

properties of the continuous-time solution. 
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1 Introduction 

The purpose of this paper. This paper is concerned with the analysis of a class of higher- 
order trajectory planning problems that are important in a wide range of contexts, ranging 
from computational anatomy to quantum control, both of which are discussed in this paper. 
The problem in its abstract formulation consists of finding an optimal curve g{t) in a Lie 
group G that generates a curve q{t) = g{t)QQ in an object manifold Q, such that q{t) passes 
near a set of fixed target points at given node times. Precise definitions of optimality and 
proximity will be provided below. 

In the higher-order variational formulation appropriate for such problems, the Euler- 
Lagrange equations split into a set of Euler-Poincare equations that hold on the open time 
intervals between the nodes, and a set of node equations that describe how to pass from one 
open interval to the next. The primary purpose of this paper is to develop a new geometric 
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understanding of the node equations in terms of conjugate momenta. Continuing from this, 
we develop a numerical algorithm for the trajectory planning problem that respects the 
geometric properties exhibited by the continuous-time solution. 

1.1 Background and problem formulation 

The problem treated in this paper fits into a classical type of problem in control theory called 
trajectory planning, or interpolation by variational curves. The task in this type of problem 
is to find an optimal curve that interpolates through a given set of points (or configurations) 
lying in some manifold, specific to the application one has in mind. Such trajectory plan- 
ning problems are relevant in numerous applications, for example, in aeronautics, robotics, 
computer-aided design, air traffic control and more recently, computational anatomy. Some 
trajectory planning applications require the optimal trajectories to possess a certain degree of 
smoothness. This requirement summons variational principles that depend on higher-order 
derivatives of the interpolation path, such as acceleration (rate of change of velocity) or jerk 
(rate of change of acceleration) etc. Properties of such higher-order variational principles 
have been widely studied, one of the earlier references being [dLR85], where an intrinsic 
formulation in terms of higher-order tangent bundles was given. [PMRRll] contains a lit- 
erature overview concerned both with mechanics and field theories, and some further recent 
developments are given in [GBHM+12, GBHRll, CdDll]. 

In some applications, for example in computational anatomy [GBHM^12] and quantum 
control [BHM12], the optimal trajectory is generated via a group action. That is, a curve g{t) 
in a Lie group G acts on a point Qq in an object manifold and generates a curve q{t) = g{t)Qo, 
that passes through a sequence of target points at prescribed times. In some instances it is 
desirable to relax the target constraints in such a way that the optimal curve does not exactly 
pass through the target points, but still passes near them at the prescribed times. This may 
be achieved by including a soft constraint in the cost functional, that is, a term penalizing 
the discrepancy between the trajectory q(t) and the targets. This leads to cost functionals 
of the following type, 

S= I l{i,t...,i^^-^^)dt + ^^Y.'^\q{t,),q.;). (1.1) 

Here, ^'^^^ are the j-th time derivatives of a curve ^(t) in the Lie algebra (the tangent space 
at the identity element of the Lie group) that integrates to the curve g{t), which in turn 
produces the trajectory in the object manifold according to q{t) = g(t)Qo. The first part of 
the cost is the integral over a Lagrangian i and is associated with the curve on the group. 
The second part sums up the squares of the distances d between the curve q(t) and the target 
points qi, at the prescribed times tj. The tolerance parameter a may be adjusted to suitably 
weight the two parts. In computational anatomy for example, a diffeomorphism group acts 
by transforming a medical image (or sub-structures thereof, such as points of interest, fibres, 
or surfaces). See [MTY02] for an overview. In this context the prescribed target images 
(or sub-structures) may not be diffeomorphically related to the initial one. That is, the 
initial and target configurations may lie on different group orbits. In general one expects 
therefore that exact matching may not be possible and works instead with a soft constraint. 
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Two of the authors have recently studied a problem of similar nature in quantum control, 
see [BHM12]. There one considers the group of unitary matrices acting on quantum state 
space, with the goal of finding the optimal experimental manipulation of the system such 
that the evolution of an initial state passes near a sequence of given target states. The cost 
functional is directly related to the required amount of change in the experimental apparatus 
over time. The introduction of a soft constraint is appropriate in this problem even though in 
this case the group action is transitive, since by increasing the tolerance parameter a optimal 
trajectories may be achieved at a lower cost. 

In a more general sense such trajectory planning problems can be thought of as inverse 
problems, where the data points qi & Q have been determined experimentally at the times tj 
and one seeks the corresponding curve g{t) in configuration space G. In this context a natural 
choice of the tolerance parameter a would be a measure of the uncertainty inherent in the 
experiment, such as standard deviation. The Lagrangian C. represents a modeling choice and 
is specific to the application one has in mind. 

This paper is concerned with trajectory planning problems of the above type. It was found 
in [GBHM+12] that the Euler-Lagrange equations split into higher-order Euler-Poincare 
equations, which hold on the open intervals between node times, and a set of node equations 
that describe the passage from one open interval to the next. These node equations describe 
the continuity properties of a set of quantities involving derivatives of various orders of the 
Lagrangian. As we will see in examples, a natural choice of Lagrangian £ leads to Riemannian 
cubics and their higher-order generalizations. This class of curves was introduced in [NHP89] 
and has since been studied in a series of papers including [CS95, CSC95, CSCOl, Noa04, 
Noa06, GGP02]. Riemannian cubics appear in a variety of applications, for example in the 
quantum control problem mentioned above, but also in computer graphics, robotics and 
spacecraft control [PR97, ZKC98, HB04, VT12]. 

Plan of the paper. In Section 2 we shall rederive the Euler-Lagrange equations for the 
higher-order variational problem by using Lagrange multipliers in a generalization of the 
symmetry reduced Hamilton-Pontryagin principle of geometric mechanics. In this approach, 
the derivation of the Euler-Lagrange equations simplifies considerably and a new geomet- 
ric interpretation of the node equations emerges. Namely, they describe the evolution of 
Legendre-Ostrogradsky momenta across the nodes, in which the highest-order momentum 
experiences a discontinuous jump related to the distance between the curve in the object 
manifold and the target points. The discontinuity can be understood in terms of a momen- 
tarily broken symmetry at the node times. However, if the object manifold is isotropic with 
respect to a subgroup action then a residual symmetry remains. By Noether's theorem, this 
residual symmetry leads to a conservation law across node times. 

In Section 3 we discuss a number of applications, including rigid body splines, macro- 
molecular configurations and quantum splines. Section 4 is concerned with the numerical 
solution of the inexact trajectory planning problem. More precisely. Section 4 describes a 
geometric discretization of the higher-order Hamilton-Pontryagin principle for inexact tra- 
jectory planning, similar to the approach given in [BRM09] for first-order systems. Our main 
motivation for the development of a geometric integrator is the exact momentum behavior 
of the discrete solution. This leads in turn to a dimensionality reduction of the search space 
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2 Geometry of the trajectory planning problem 

We start with the statement of the problem considered here. One aims at steering from an 
initial point Qo in some object manifold Q along an optimal trajectory q{t) that evolves via 
the action of a Lie group G. That is, q{t) = g{t)Qo, where the right hand side denotes the 
action of g{t) on Qq and the curve q{t) lies in the G-orbit of Qo- 

The optimality condition is given in terms of a function £ : fcg — )• R defined on the fc-fold 
Cartesian product kg of the Lie algebra g, which measures the cost of the transformation 
g{t), and a distance function d : Q x Q ^ H. As we shall see, the integer k determines 
the degree of smoothness of solution curves. The optimal curve q{t) is required to pass near 
prescribed target points Qt. at prescribed node times ti for i = 1, . . . , /. This is formalized 
by including a squared distance term d'^{g{ti)QQ, Qt-) in the cost functional, for each i. Thus, 
the cost functional S* : C(0) — > R, where C{g) a suitable space of 0- valued curves (see below) 
is defined by 

S[i\ ■■= f ^(e, • • • , e^'-'^) di^ + ^ E d\g{U)Q,, Qu). (2.1) 

The notation ^^^^ is shorthand for d^^/dP , the quantity cr is a tolerance parameter, and the 
curve g(t) originates at the identity g{0) = e and satisfies g = -^\ _r,exp(£^)5f. We write Rg 
for multiplication by g from the right and TRg for its differential, thus g = TRg^. Variations 
are considered in the space C{q) consisting of curves ^(t) : [0,tj] — )■ Q whose restrictions to 
open intervals (tj_i, tj) for i = 1, . . . , / are C"^^'"^ and whose j-th derivatives ^^^^ are continuous 
on [0, ti] for j = 0, . . . , A; — 2. We also assume that initial values of ^'■■'''(0) for j = 0, . . . , /c — 2, 
are given. 

This type of trajectory planning problem is familiar, for example, from image registration 
in computational anatomy, where one typically thinks of Qq as a template shape being 
deformed by a curve of diffeomorphisms g{t), in turn generated by the time-dependent vector 
field ^(t) [YoulO]. At times tj the resulting curve in shape space passes near the given target 
shapes Qt^, the parameter a determining the proximity of the passage. In this case, the 
Lie group G of diffeomorphisms is infinite dimensional. However, in the present paper we 
will restrict ourselves to the case of finite-dimensional Lie groups and object manifolds. 
A natural finite-dimensional instance for illustrating these ideas arises in quantum control 
[BHM12], where quantum state vectors evolve under the action of the unitary group. The 
generator curve ^(t) in this case corresponds to the Hamiltonian operator, which is controlled 
in experiments. 

2.1 Euler Lagrange equations via Lagrange multipliers 

The Euler-Lagrange equations characterize solutions to Hamilton's principle 5S = Q and 
were derived in [GBHM+12]. The equations split into a set of higher-order Euler-Poincare 
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equations on the open time intervals between the node times and a number of node equa- 
tions describing how the solution evolves across the nodes. In previous formulations of this 
problem, one must find the variation 6g{ti) that is produced by a (time dependent) variation 
SC,{t). This can be facilitated by taking advantage of Lagrange multipliers in an equivalent 
variational formulation that we describe now. As this paper demonstrates, the new approach 
also provides a geometric interpretation of the node equations and furthermore suggests a 
geometric numerical procedure for the solution of the problem , see Section 4 below. 

The method of Lagrange multipliers involves enlarging the space on which the dynamics 
happen. We define the cost functional 5* on some space of curves C{G x kg x kg*), 



s[9,e,.--,e-' ■■' 



,/^ 



,/^ 



fc-ii 



1 I 
+ J2^f^'^r'-n dt+—j2d'i9iu)Qo,Qu). (2.2) 



fc-i 



r=l 



j=l 



Again some technical assumptions about the space of curves are needed. Namely the curves 
are C^ when restricted to the open intervals (tj_i,tj) for i = 1,...,/, and (7,^°, . . . ,^'^~^ 
are continuous on [0,t;]. We also assume g{0) = e and given initial values ^■'(O) = C,o, for 
j=0,...,fc-2. 

Before we take variations of S it is useful to introduce the cotangent lift momentum map 
J^ : T*Q — )■ 0* associated with the action of G on Q (see []\1R03], Chapter 11, for more 
information). This map is defined to satisfy 



(a 



qi 



iQ{q)) = (^'^K),0, for any a, G T*Q, ^Eq. 



(2.3) 



where we have used the notation ^^(g) := ^I^^q^^^^ ^^'^ (•'•) for the respective duality 
pairings. We also introduce the shorthand di(i(gi, 52) G T*Q to denote the exterior derivative 
of the distance function d with respect to the first entry, and di(i(tj) := did{g{ti)QQ,Qt.). 
Integrating by parts and using (2.3) we obtain 



6S 
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k-2 



(2.4) 



r=0 
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where we set rj := TRg-iSg and defined fi^{t~) := limift^ f^^it) as well as /U^(t^) := lim^t^ f^^{t)- 
The last line above could have been omitted since by assumption 77(0) = and 6^,-^ = for 
j=0,...,fc-2. 

We can now read off the Euler-Lagrange equations. On the one hand, for t in any of the 
open intervals (tj, tj+i), i = 0,...,/ — 1, we have 

/i° + ad^o /i° = 0, (2.5) 

TRg-ig - e° = 0, (2.6) 

C-'-C = 0, (r = l,...,fc-l) (2.7) 

f,r + ^r-i__p_^Q^ (r = l,...,A:-l) (2.8) 

/-^ - ^ = 0- (2-9) 

On the other hand, the node equations are given by 

/(t;) - /i°(t+) + ^ j«(dirf(t,)) = 0, (s = 1, ...,/- 1) (2.10) 

/i'^(t7) - fi'itt), (r = 1, . . . , fc - 1; s = 1, . . . , Z - 1) (2.11) 

/(to + ^ J^(dirf(tO) = 0, (2.12) 



a 



fi'iti) = 0. (r = l,...,fc-l) (2.13) 

Remark 2.1. There are 4 versions of the action functional, which are all relevant in appli- 
cations. The one above can be called the left-action, right-reduction version since g{t) acts 
on Qq from the left, while C,^ is the right-reduced velocity C,^ = TRg-ig (see Section 2.2 below, 
for more details on reduced variables). There are the following three other cases. 

(1) The right-action, right-reduction case with action functional 






i{e,e,■■■,e-') + {^^',TRg-.g-e) 



r=l -' i=l 

The variation of the penalty term in this case changes according to 

^Sd\g{ti)-'Qo,Qu) = -d{t,) (Ad;(,,)-i J{d,d{t,)),r]{t,)) 
This means that (2.10) and (2.12) are replaced by 

Atj) - Att) - ^ Ad;(,^)_. j^idAts)) = 0, (. = 1, . . 
Ati)-^Adi^t,)~^j^{d,d{t,)) = o. 



i-i) 
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(2) In the right-action, left-reduction case, the action functional is 



S[G,E',...,E'-\m',...,m'-'] 



l(E^E\ . . . ,E'-') + {m^TLc-.G -E') 



fe-l n -, / 



r=l 



+ X^K, s-i - s^) dt + ^Yl d'iGitr'Qo, Qu 



2a2 



where we wrote Lq for multiplication by G from the left and TLq for its differential. 
However, this is equivalent to the left-action, right-reduction case by identifying 

G = g-\ E^ = -^^,...,E''-^ = -^''-\ m° = -/,..., m'^-i = -/-I (2.14) 

and setting i = I o n, where k, : kg —^ kg is multiplication by —1. 

(3) By the same token the left-action, left-reduction case with action functional 



5[G,S°,...,S^-\m°,...,m'^-i] 

'0 

fc-i 



ti 



l{E\E\...,E''^) + {m\TLG-^G~E') 



fc-i -| / 

+ J]K,S^-i-S'-) dt+—Y,d\G{t;)Q,,Q,:). 

r=l -' i=l 

can be mapped to the right-action, right-reduction case. 

In the analysis that follows we will largely restrict ourselves to the left-action, right-reduction 
case. Anything we say can be transferred to the remaining three cases by applying the modi- 
fications listed above. 

2.2 Euler Poincare equations 

From (2.5)-(2.9) it follows that on open intervals (tj, tj+i), i = 0, ...,/ — 1, 

fc-i 



d 
~dt 



/ j=o '' 



(2.15) 



This is a k-th order Euler-Poincare equation for a system that exhibits right-invar iance. This 
type of equation was derived in [GBH]\I+r2] from a variational perspective and in [CdDll] 
from the Hamiltonian one. We now explain its appearance in the inexact trajectory planning 
problem from the viewpoint of invariant Lagrangians, starting with some necessary definitions 
that can be found in [CMROl] and [GBHRll]. 

The k-th order tangent bundle Tq : T^^^G — )■ G is defined as a set of equivalence classes 
of curves as follows: Two curves gi{t) G G, i = 1,2, are equivalent, if and only if their time 
derivatives at t = up to order k coincide in any local chart. That is, f?! (0) = (72 (0), 
for Q < j < k. The equivalence class of a curve g{t) is denoted by [(?] (g)' *-'^ formally as 
(5^(0), 5'(0), . . . , f7*^^^(0)). The set of all equivalence classes of curves emanating from Qq E G is 
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written as Tg^ G. Then T^'^^G := IJoeG-^3 ^ ^^ ^ fibre bundle over G with projection map 
Tq : [g]„rQ\ ^-^■ g{0). Note that a curve ^'(t) defines, at each time t in its domain, an element 

[9]% ■= [h]hl) by setting h{s) := g{t + s). 

It is convenient to represent T^'^^G using the trivialization map that makes use of the 
right group multiplication (analogous constructions exist using the left multiplication map). 
Let g{t) be a representative of [g]„(Q\ and define ^{t) := TRg~ig{t) . The trivialization map 
T^^-'^G — )■ G X /eg is given by 

The reduction map (3k : T^'^'^G — > fcg is obtained by omitting the first entry. 

It is a well known (see, for example [MR03], Chapter 13) that the first order [k = 1) Euler- 
Poincare equation appears when the Euler-Lagrange equation for a Lagrangian TG — )■ R 
with symmetry is written in terms of the reduced velocity vector ^{t) = TRg-i(^i^g{t) . The 
higher-order Euler-Poincare equation appears in a similar fashion when computing the Euler- 
Lagrange equations for an invariant k-th. order Lagrangian L : T^^^G — )■ R. More precisely, L 
is called right-invariant if L{\yg] L^ = L(^[gh] ,Lf^y The definition of left-invariance follows 
analogously. Let L be a right-invariant Lagrangian and consider Hamilton's principle, SJ^ = 
0, for 

J[g]= [ L{g{t),g{t),---,9^'\t))dt, (2.16) 



where variations are taken with respect to fixed end points up to order k — 1, that is, 
5g^^\a) = 5g^^\b) = 0, for j = 0, . . . fc — 1, in any local chart. The Euler-Lagrange equations 
can be written in terms of the right-reduced velocity vector, which leads to the k-th order 
Euler-Poincare equations [GBHM+12] 

/ j=0 '' 

with reduced Lagrangian £ : A;g — )■ R, 

It is now straightforward to see from the viewpoint of invariant Lagrangians why the Euler- 
Poincare equation (2.15) must characterize optimal curves in the trajectory planning problem 
on open time intervals. Indeed, fix < i < / and suppose ^(t) is a local extremum of (2.1) with 
integral curve g{t). H C,\{ti,ti+i) is not a solution to the k-th order Euler-Poincare equation one 
can find a variation 5'e(t) keeping fixed g(ti) and g(ti^i), such that S /^.*^^ L(g, g, . . . , g^^^) dt 7^ 
with Lagrangian L := io p^. By consequence, ^£{t) := TR -i,^^gs{t) is a variation of ^ such 
that 6S y^ ioT 5* as in (2.1), a contradiction. 
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2.3 Geometry of multipliers 

Furthermore we observe from (2.5)-(2.9) that on open intervals (tj,tj_|_i), i = 0, . . . ,1 — 1, 

"'-E'i-D'Ij;,^ (r = 0.....t-l). (2,18) 

In order to discuss the geometric meaning of these identities, we recall from above that 
the trajectory planning problem (2.2) on open intervals reduces to a problem of the type 
6J^ = with J' of the form (2.16). In particular, equations (2.5)-(2.9) and therefore (2.18) 
are obtained by taking suitably constrained variations of 

fc-1 



+ J2(f^\r'-n 



r=l 



dt (2.19) 



This variational principle is a higher-order generalization of the reduced Hamilton-Pontryagin 
principle of first order mechanics. In first order mechanics, this principle provides a unified 
treatment of the Lagrangian and Hamiltonian descriptions of invariant mechanical systems 
on Lie groups (see [YM06] for a detailed discussion^). In particular, the Legendre transform 
connecting the two descriptions is revealed by the variational calculus. This remains true 
for higher-order mechanics. Indeed, (2.18) can be recognized to be the reduced Legendre 
transform that appears in [GBHM+12, GBHRll]. While we found (2.18) from a variational 
approach, these references take as starting point [dLR85], where a coordinate free descrip- 
tion of the higher-order Legendre transform on manifolds was given. We briefly review this 
approach here. 

The Legendre transform of higher-order mechanics, given in [dLR85], is a map Leg : 
2"(2fc-i)(^ _v. T* {T^^~^^)G . If Leg is a diffeomorphism (that is, L is hyperregular) it connects 
the Lagrangian and Hamiltonian descriptions just as in the first order case. With respect to 
the right-trivializations 

- 1)0* (2.20) 



rj.{2k~l)Q ^Gx{2k- 


-2)0, T*(T('=-^)) = Gx (A:-2)0X ( 


it is given as [GBHRll 




Leg:((7,e°,. 


..,e'-')-^i9,e,...,e-\f^',...,f^'-' 


where 


f^'=~i2i ^yfJl {r=o,...k 

"^ ^ 'dp S^-'+i ^ 

3=0 ^ 



I). 



The same equations were seen in (2.18) to emerge from the Hamilton-Pontryagin principle 
(2.19). This means that, as for first order, the higher-order Hamilton-Pontryagin principle 



^ There is a close connection between the Hamilton-Pontryagin principle and Dirac structures, an aspect 
we do not enter into in the present paper. See [YMOG] for more information. 
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contains both Lagrangian and Hamiltonian descriptions of higher-order mechanics and pro- 
vides a unified framework for both views. To obtain the Lagrangian description one may 
ehminate /i°, . . . ,/i'^~^ from (2.5)-(2.9) using (2.18). The resulting equations are the trivial- 
ized fiow equations of the Lagrangian vector field, an element of X(T*^^'^~^^G) [dLR85]. On 
the other hand, if (2.9) can be solved for C,^~^ (this is the case, for example, when L is hy- 
perregular) then (2.5)-(2.8) are the trivialized fiow equations of the Hamiltonian vector field 
Xh e X(T*(T('=-^)), which solves [(1LR85] 

ix„u = dH, (2.21) 

where u is the canonical symplectic form on T*{T^^^^'>G) and H : T*(T^''^^^G) — > R is given 
as 

fc-i 

H{g, e,..., e-\ /i°, • • • /-') = E ('"'■' ^^) - ^(^°' • • • ' ^'"') (2.22) 

r=0 

with respect to the trivialization (2.20). By consequence of (2.21) the fiow map Ft : 
T* {T'^^~^^) — )■ T*(T*^'^~^)) of the Hamiltonian vector field preserves the symplectic form u. 

For later reference we point out how this can be seen alternatively from the Hamilton- 
Pontryagin principle. This is a generalization to higher order of a standard argument (see 
for example [BRM09], Section 3, for the first order case). If we omit end point constraints 
on the variations oi J in (2.19), the integration by parts contributes boundary terms to 5 J 
(cf. (2.4)), 

h k—2 h 

6J= I ■■■ dt+ {^i\TR,-.6g)\l + Y.{^'''^'^^^'')t= I ■■■^t+0{6x)t, (2.23) 

Ja j,_Q Ja 

where 6 in the second equality is the canonical one-form on T*{T^^'^^G) and we defined 
5x{t) to be the curve in TT*{T^''~^^G) whose trivialization corresponds to the variations 
{TRg-i5g, 5^^, . . . , 5^^''~'^\ SfjP, . . . , (J/x^^"^-'). If we restrict the variations to solution curves 
of (2.5)-(2.9), we may just as well express ^7 as a function of initial conditions i/jnitial " 
T*(T^^^^^G) — !■ R. The integral part of (2.23) then vanishes and we obtain 

5J = dj;^iti^^{5x{a)) = {Fl.e - e){5x{a)). 

Therefore, F^_J) = 6, and taking an exterior derivative yields the desired identity F^_^uj = u. 

2.4 Momentum conservation and Noether's theorem 

Another important point to notice from (2.5)-(2.9) is that Ad*yU° is a conserved quantity. 
Here Ad* is the map dual to Ad : G x g — )■ g given by {g, C,) i— )■ Ad^ C, := TLgTRg-iC,. Indeed 
by (2.5) 

^ Ad; /i° = Adlifi' + ad^o /) = 0. (2.24) 
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In the context of first order Euler-Poincare equations a similar momentum conservation is 
due to the invariance of the Lagrangian with respect to group multiplication operations. This 
is an instance of Noether's theorem, which roughly speaking guarantees that the momentum 
map associated with the action of a symmetry group is preserved (see for example [MR03], 
Chapter 11). We now show that the situation is similar for (2.24). The right action R oi G 
on itself, 

R:GxG, {h,g)^Rg{h) = hg, 

can be lifted to an action on T^'^^^^G, 

y(.-i)^ : T(-^)G X G -, T^'-'^G, {[h]%l\g) ^ T^'-'^R, {[h]i%l^) = [hg][\-l^. 

This action can be lifted to its so-called cotangent lifted action ([MR03], Section 12.1) 

rp*rp{k-i)j^ . ^^(^(fc-i)^) X G ^ T*(T('=-i)G), 

given in trivialized form as 

T*T('=-i)i?,(/^, e,... e'=-^ /, . . . /-^) = {hg~\e, .., e-', /, . . . /-^). 

It is apparent that the Hamiltonian (2.22) is symmetric with respect to this group action. 
By Noether's theorem the associated momentum map is conserved. 

What is this momentum map? By appealing to standard formulas ([MR03], Section 12.1) 
we find that the momentum map J : T*{T^'^^^^G) — )■ Q* is 



J(^,e°,...e'-',/,.../-')=Ad>' 



with respect to the trivialization (2.20). The conservation law observed in (2.24) thus arises 
from the right-invariance of the Lagrangian (respectively, the Hamiltonian) via Noether's 
theorem. 

The conservation law can also be obtained from a variational perspective. This is well 
known in first order mechanics, and it is also the case in higher-order mechanics. We take 
a solution of (2.5)-(2.9) on the time interval [a, 5] and vary it according to 5g = TLgV for 
z/ e 0. For J as in (2.19) we have (cf. (2.23)) 

5J = Q= (/, TRg^.TLgu) I = (Ad; /, u)\l. (2.25) 

The same argument holds after replacing the upper boundary b by any b' G [a,b]. Since u 
was arbitrary we conclude that Ad^/x^ is conserved along a solution of (2.5)-(2.9). 

2.5 Node equations 

The remarks above concerned equations (2.5)-(2.9) on the open time intervals between nodes. 
We now come to the node equations (2.10)-(2.13). These specify the evolution across node 
times of the Lagrange multipliers /i*", which we interpreted above as the reduced Legendre 
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momenta of the system. More specifically, the momenta fi"^, r = 1, . . . ,k — 1 are continuous 
on [0,ti\, while the 0-th momentum /i° experiences jump discontinuities at the nodes. If the 
Lagrangian i is hyperregular we can conclude that g e C^'^~^([0,ti]), that is, g is {2k — 2) 
times continuously differentiable on [0,t;]. Furthermore, the node equations specify terminal 
values for the curves /i*", r = 1, . . . ,k — 1. 

For a, 6 G R define la<b to be equal to 1 if a < 6 and otherwise. We can now prove the 
following theorem. 

Theorem 2.2. For t in any of the open time intervals (t^j^s+i) as well as for t G {0,ti}, 

/(t) = -^Ad;(,)-. f ^l,<,^rf(t,)Ad;^ J«(did(t,)) j . (2.26) 

Proof. At final time t = ti (2.26) clearly holds because of (2.12). Since Ad*/i° is conserved 
on open intervals it follows that for t G (^5,^^+1)5 

/(t)=Ad;(,)-.Ad;(,^^^)/(t;_,i). 

We can now obtain (2.26) by induction over the open time intervals, noting that at each node 
t = ta a. term — ^J'^(di(i(ts)) gets added on. ■ 

In order to formulate the following corollary we use the notation Qq for any point q & Q 
to denote the Lie algebra of the isotropy subgroup of that point, G q := \^g E G\gq = q^ . In 
particular pqi^q) = for any p E Qq. 

Corollary 2.3. For a solution of (2.5) -(2. 13) we have, fort in any of the open time intervals 
[ts, ts+i) as well as for t G {0, ti}, 

(/(t),p) = for all p e Qqt^t)- (2.27) 

In Section 4 we will develop a geometric algorithm that inherits an exact version of this 
Corollary. This implies that the numerical search for the optimal initial value of /i° can be 
restricted to the subspace of q* that annihilates Qq^. 

Proof. For t and p as in the statement of the corollary it follows from Theorem 2.2 that 

I 



'yit),p) = -^Yl i*<*= (c^(t.) j«(dirf(t,)), Adg,^ Adg^t)-^ p) 

s=l 

= ^5^ U<u {d{ts) dirf(t,), (Ad^,^ Adg(i)-i p)q {q{ts))) = 0, 



where we used (2.3) for the second equality and noted that Ad^^^^ Adg(t)-i p G 5q{ts) for the 
third. ■ 
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2.6 Residual conservation law after partial symmetry breaking 

A physically intuitive perspective on Corollary 2.3 is to understand it as a residual conser- 
vation law after partial symmetry breaking. We can see the sum in (2.2) as the integral over 
a time-dependent potential function V : [0,ti] x G ^^ M, given by 

1 ' 
^(^' 9) = 7^Y.^^^- ^^) d\9it)Qo, Qu) (2.28) 



i=l 



This potential produces instantaneous singular forces at node times tj which impart the jump 



discontinuities on the otherwise conserved momentum J = Ad* jj' 



* ,,0 



J{tt) = Jit;) + ^ Ad;(,^) J'^idAts)). (2.29) 

This is because the presence of this potential breaks the G-invariance of the variational 
problem, however a residual symmetry remains. Clearly, multiplication of g from the right 
by an element h G Gq^ leaves V invariant. An adaptation of the argument surrounding 
equation (2.25), restricting u to the subspace 0q„ C q, then leads to 

o = (Ad>°,z.>i;. 

for any t E [0,ti]. Moreover, (2.12) guarantees that (Ad*^^ ^ /i'^(t;), z/) = 0. Therefore, 
(Ad*((-) /i°(t), z/) = for any t G [0,t/] and u G Qq^, which is equivalent to Corollary 2.3. 

3 Applications 

In this section we discuss a number of examples that summon the inexact trajectory planning 
problem. We start with a discussion of Riemannian cubics on general Riemannian manifolds 
and specifically on Lie groups with invariant metrics, since this type of curve appears for 
a certain natural choice of Lagrangian. We then treat the rigid body, finite-dimensional 
quantum systems and molecular strands. 

3.1 Riemannian cubics 

Let (M, 7) be a Riemannian manifold with metric 7, whose norm we denote by ||.|L. The 
notion of straight lines in Euclidean space generalizes to geodesies in M. These are curves 
x{t) G M that satisfy DtX = 0. Here we introduce the notation DfX = V^i;, where V is the 
Levi-Civita connection for 7. In local coordinates the geodesic equation is given by 



A;fc 



X 



+ T'^^x'x' = 0, 



«j 



where F^- are the Christoffel symbols of the Levi-Civita connection. We define the vector 
bundle isomorphism over the identity, b : TM — )■ T*M, which maps v G T^M to v^ = 
l{x){vx, ■)• The inverse map is denoted by jj. 
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The physical interpretation of straight lines in Euclidean space follows from Newton's 
second law x = F. Here F is an external force acting on a particle of unit mass that follows 
the trajectory x{t). By analogy one can often understand {DtxY as a generalized force acting 
on a physical system whose time evolution is represented by a curve x{t) in configuration 
space M. Consider the Lagrangian L : T^'^^M — )■ R, 

1 2 

L{x,x,x) = - ||Aa;|L, (3.1) 



2 



whose corresponding action functional S = J Ldt measures the square of the L^-norm of 
the external force. Hamilton's principle 6S = for variations with fixed boundary velocities 
x{a) and x{b) leads to the Euler-Lagrange equation [NHP89] 

Dfx + R{DtX, x)x = 0, 

where R is the curvature tensor, R{X, Y)Z := VxVy2' — Vy^ xZ — 'S/[x,y]Z for any vector 
fields X,Y,Z G X(M). Solutions to this equation are called Riemannian cubics and were 
introduced in [NHP89]. 

When the manifold is a Lie group M = G we call a Riemannian metric right (or left) 
invariant ii j{g){vg,Wg) = j{gh)(TRhVg,TRhWg), or 'y{g){vg,Wg) = 'y{hg)(TLhVg,TLhWg) 
respectively, for all g, h & G and Vg, Wg G TgG. If the metric is right (or left) invariant, then 
the Lagrangian (3.1) can be reduced to a function £ : 2g — )■ R [GBHM+12] 

iie,e) = l e'±adU° ', (3.2) 

Z 7 

where we introduced the operation adlp = {a(Kp")K This can be seen from the fact (e.g., 
[GBHM+12]) that for a curve g{t), whose right (or left) reduced velocity is given by ^(t) G Q, 

A^= (e + ad|e)^, or Dtg = g(i-a.dl^y 

In the following we will discuss a number of physical systems whose configuration spaces 
are Lie groups and whose equation of motion in the absence of external forces is given by 
Dtg = for some right (or left) invariant Riemannian metric. The Lagrangian (3.2) is then 
one natural choice for i in the inexact trajectory planning problem since optimal curves 
minimize (in the L^ sense) the amount of external forcing necessary to achieve them. It is 
clear from equations of motion (2.5)-(2.11) that the solution g{t) is a Riemannian cubic on 
open intervals, and twice continuously differentiable on the whole time interval [0,t;]. Such 
curves are called Riemannian cubic splines. 

Remark 3.1. Without going into the mathematical details we point to a probabilistic inter- 
pretation of Riemannian cubics. See also [VT12J, where a closely related idea is discussed 
in the context of stochastic modeling of biological growth. Let G be a Lie group with right- 
invariant metric 7 and let Ci, i = 1, . . . ,d be an orthonormal basis of the d-dimensional Lie 
algebra q. Consider a curve g{t) G G, whose right-reduced velocity ^ = TRg-ig satisfies the 
following Ro stochastic differential equation 

d 
d^ = -adJ^dt + avr^diyiCi, (3.3) 
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where Wi, i = 1, . . . ,d, are independent Brownian motions (see, for example, [0ksO3j Section 
2.2, for a definition) and aw ^ ^- Suppose the (noisy) data is given in a vector space V 
equipped with an inner product, whose norm we denote by \\-\\v- The noise distribution 
is assumed to be Gaussian, that is, the probability density function has the form p{Q) ~ 
exp(— ^IIQ— Ollv'); whereQ is the true state of the system andan G R. Suppose experiments 
at times ti, % = 1,...,/ measuring the trajectory g{t)Qo have given results Qt,.. Then the 
minimization of 



S= ri{^,Odt+^J2\\9it^)Qo-Qu\\l 

Jo '^^n ._i 



with I as in (3.2), can formally be understood as the maximization of the (logarithm of the) 
probability of the path git), given the measurements. Alternative models of stochastic forcing 
will typically lead to minimization problems of the same type, but with different choices of i 
(see Remark 3.2 below). 

3.2 Rigid body splines 

Let the Lie group G be the set of rigid rotations 5*0(3), and let Q be the unit sphere S'^ G R^. 
We use vector notation for the Lie algebra 50 (3) of the Lie group of rotations 5*0(3), as well 
as for its dual 5o(3)*. One identifies so(3) with R'^ via the familiar isomorphism 



: R^ ^ so(3), ri=6^fi:=ri= a -c, (3.4) 





called the hat map. This is a Lie algebra isomorphism when the vector cross product x 
is used as the Lie bracket operation on R^. The identification of 5o(3) with R^ induces an 
isomorphism of the dual spaces so(3)* = (R^)* = R'^, with the dot product as duality pairing. 
Let 7 be a left-invariant Riemannian metric on 5*0(3). This defines an inner product on so (3) 
which can be expressed as 

7e(f2l,f22) = rii-m2 

for a symmetric, positive definite matrix I. The geodesic equation Dtg = in Euler-Poincare 
form is 

i7 + ri(ri X m) = 0, g = g^. (3.5) 

Consequently, the Lagrangian (3.2) takes the form 

i{n°, n^) = -{n^ + r^n'^ x m")) ■ 1(^2^ + r^(f2° x if2°)) 

= l\\n' + r\n'xmm^^^. (3.6) 
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Consider the inexact trajectory planning problem in the left-action, left-reduction form. That 
is, suppose an initial point xq and targets Xj^, . . . x^^ G S"^ are given, as well as a tolerance 
parameter a. We seek the minimizer of 

s[g,n',n\^^^^'] = T ^ Hfi^ + ri(f2° x m°)||l(3) + {f^\g-'g - n') 



1 

+ (^1,17° - f^i) + — 5^ \\g{t,)^o - ^u 

i=\ 



(3.7) 



The physical interpretation is as follows. The group of rigid rotations, 5*0(3), is the config- 
uration manifold of a rigid body constrained to rotate around a fixed point. In the absence 
of external torques the motion is governed by the geodesic equation (3.5), where I is the 
moment of inertia tensor (see, for example, [HoUl], Section 2.4). The resulting curve g(i) 
describes the orientation of the rigid body relative to a space-fixed reference frame. 

Suppose the motion of a rigid body (with or without external torque) is partially observed 
in an experiment. At discrete times tj the direction of a particular body fixed axis is measured 
in the space-fixed frame, generating a sequence of outcomes x^^ G S*^. Therefore, if xq is the 
initial direction of the axis and git^ describes the rigid body motion, then git^'x.Q — x^. = 0, 
up to measurement error. One would like to model the trajectory 5'(t), taking into account 
this information. The action functional (3.7) encodes one such model, yielding the curve 
g[t) of minimal external torque (in the 1? sense) that is consistent with the experiment. A 
natural choice for the parameter o"^ is to set it equal to the variance of the measurement. 
Example simulations are shown in Figures 3.1 and 3.2. These were generated by the numerical 
algorithm discussed in Section 4. 




Fig. 3.1: Trajectory planning on the sphere through 
rigid rotations. The moment of inertia tensor was 
taken to be the identity matrix and the tolerance pa- 
rameter was set to cr = 0.025. The given data points 
are represented as black dots with uniform time sep- 
aration. The curve shown is g{t)x.Q for the optimal 
curve g{t), generated using the algorithm developed 
in Section 4. The coloring represents the local speed 



along the curve in SO (3), that is, 
large, white is small). 



\n'> 



lso(3) 



(red is 



If the observer is instead moving with a body fixed frame and measuring a space fixed 
direction, then (7(tj)~^Xo — Xj. = 0, up to measurement error. In this case the inexact trajec- 
tory planning problem presents itself in the right-action, left-reduction form. Evidently, the 
formalism presented in this paper applies to any (sufficiently smooth) choice of Lagrangian. 
For example. 



[ri°, n^) = -f2° ■ {m^ + ri° x m°) 
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Fig. 3.2: Optimal curve g(t) in the group 5*0(3) of rigid rotations corresponding to the data points shown in 
Figure 3.1 and generated using the algorithm discussed in Section 4. The moment of inertia tensor was taken 
to be the identity matrix and the tolerance parameter was set to a = 0.025. In this figure a given element of 
5*0(3) corresponds to a point on the radial line along the rotation axis, at a distance from the origin equal to 
the rotation angle. The radius of the inner sphere is tt and the radius of the outer sphere is 2tt. The center, 
marked by a cross, and the boundary of the outer sphere thus both represent the identity matrix. 



leads to optimal curves g{t) with minimal work done by external torques. 

Remark 3.2. We mentioned in Remark 3.1 that the minimization of (3.7) is related to a 
certain inverse problem given the stochastic evolution (3.3). Alternative stochastic models 
lead to forms of C. different from (3.6). For example, 

dn = -r\n x in)dt + av^dw, 

where dW = (dW^"*^, dVT^, dW^^)"^ is a vector of independent Brownian motions, leads to 



with \\.\\ being the Euclidean norm. 



\n^ + r\Q.° xin^ 



3.3 Quantum splines 

In this section we give an overview of a quantum mechanical application presented in [BHM12], 
where more details can be found. We consider an n + 1 level quantum system with Hilbert 
space Ti = C""*"^. Quantum state space is the n-dimensional complex projective space CP", 
which is defined as C""*"^ — {0} modulo the equivalence relation \ijj) ~ X\ip) for any complex 
number A 7^ 0. The standard geodesic distance is given by 



d{ip, 0) = 2 arccos 



(^10) (01^) 
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Here we used Dirac notation to write lip) for a vector in C"'^^ and {ipl for its Herniitian 
conjugate. The evolution of a quantum state is given by the Schrodinger equation 

^|V^) = -ii^|^), 

where the Hamiltonian operator iJ is a Hermitian matrix. Equivalently one may express this 
as a differential equation on the unitary group U{n + 1) by 

U = -iHU, |^(t)) = t/(t)|^(0)). (3.8) 

One may assume without loss of generality that the anti-Hermitian matrix iH is of zero trace, 
therefore U{t) is in the special unitary group SU{n + 1). This is due to the fact that the 
trace contributes a complex phase factor to the state evolution, which can be neglected in 
projective terms. Notice that —iH is skew-Hermitian and trace free and therefore lies in the 
Lie algebra su{n + 1). 

Suppose a time dependent Hamiltonian H{t) is controlled in an experiment whose goal 
it is to steer a system from initial state \iIjq) through states lipt-) at times tj, for j = 1, . . . ,1. 
One would like to achieve this trajectory with least change to the experimental apparatus. 
We formalize this requirement using an inner product on su(n + 1), 

7e(A,S) = -2tr(Afi). (3.9) 

The cost we associate with a time dependent Hamiltonian H{t) is J ^'~ff.{H,H) dt. We note 
that 7e extends to a bi-invariant Riemannian metric 7 on SU{n+l), which means in particular 
that ad' = — ad. The Lagrangian i{H^, H^) = ^"feiH^, H^) is therefore equal to the reduced 
Lagrangian (3.2) for the metric 7 on SU{n+l). The total cost functional takes the left-action, 
right-reduction form 

S[U,iH^:iH\M^,M^] = r ti{H^H^) + /m^,UU-^ - iH^\ + /M\iH^ ~ iH^\ dt 

I 



^J2d'iUit,)^Po,^Pu 



2a2 



=1 



The tolerance parameter a can be used to trade off quality of matching against a reduced 
cost associated with change in the Hamiltonian over time. For more details, see [BHM12]. 

3.4 Macromolecular configurations 

Equilibrium configurations of macromolecular structures and of DNA in particular can be 
modelled using the classical theory of elastic rods. See [SH94, KC06] for examples of this 
approach. In this section we formulate an inexact trajectory planning problem in this context. 

We start by describing how the configuration of an elastic rod can be described by a 
position curve r(s) in R^ and a curve R{s) in the group of rigid rotations, 5*0(3). Here, 
s G [0, 1] parameterizes the cross sections of the rod along its length, whereby for a given 
value of s the vector r(s) points to the center of mass of the respective cross section, as 
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seen in the lab frame. Let ej(s), i = 1, 2, 3 be an orthonormal frame such that ei(s) and 
62(5) point along the principal axes of the moment of inertia tensor of the cross section. We 
will refer to this frame as the body-fixed frame. R{s) is the rotation that transforms the 
initial frame at s = (which we assume to coincide with the lab frame) to the body-fixed 
frame at s. Therefore the configuration of a macromolecule can be described by a curve 
g{s) = {R{s),r{s)) in the special Euclidean group, SE{3), originating at the identity. The 
group multiplication rule of SE{3) is 

{Ri,ri){R2,T2) = iRiR2,Rir2 + ri). 

The Lie algebra 5e(3) consists of elements (fi, v) with Q G 5o(3) and v G R.^. Applying the 
inverse of the hat map (3.4) to fl we can represent 5e(3) as R^. The ad operation becomes 

ad^i ^2 = ad(nj,vi)(f^2, V2) = {^i x fig, fii x V2 - fia x Vi). 
If we identify the dual 5c(3)* with R^ using the standard dot product as duality pairing, then 

ad(n^v)(A*) ^) = (~^ X /^ — V X a, —ft x a). (3.10) 

The body-fixed velocity vector pertaining to a configuration {R{s),r{s)) is defined as 

^ = g-'g = iR-'R,R-H) = {n,R-'r) 

where a superscript ■ means derivation with respect to s. In vector notation we can set 
V = R-^t so that $, = {n, v) G R^. 

A macromolecule that is experimentally constrained to assume a configuration with given 
final rotation and displacement g{l) = (i?(l),r(l)) will relax into an equilibrium state that 
minimizes potential energy with respect to all possible configurations respecting the con- 
straint. For the case of DNA the authors of [KG06] propose to model this effect using the 
Lagrangian L : TSE{3) — ?■ R whose left-reduced form / : se(3) — ?■ R is given by 

The 6x6 matrix K is symmetric and positive definite and encodes the various stiffness 
properties. The double helix structure of DNA means that the equilibrium configuration for 
unconstrained end points retains a number n of rotations along its length. This is expressed 
by the vector 

^7^^ ) , (3.11) 

where we assume without loss of generality that the equilibrium configuration for uncon- 
strained end points is oriented along the spatial z-axis and has unit length. 

In the case of constrained final rotation and displacement g{l) = {R{l),r{l)) the equi- 
librium configuration minimizes the action functional S = J^ /(^) ds. Hamilton's principle 
6S = leads to Euler-Poincare equations 

^ = K-^a.dlK{^-z) = adl{^-z), (3.12) 
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with ad* operation given by (3.10) and the operation ad' defined by adl ^2 •= ^~^ ^^l ^$,2- 
This equation of motion can be used to design second order Lagrangians to model non- 
equilibrium states of the DNA. For example, let us set 

e(e,e) = i(«' - adt.(«° - z)) . A-(«' - ad^,.(«° - z)) 



1 - ....o^, ^ 

K 



^'-ad^o(^o-z) 



Suppose an experiment measures the position of the center of mass r{si) at a number of 
parameter values Sj (i = 1, . . . , /) as well as a body fixed direction, say 63(54). The space of 
measurement outcomes is Q = 5*^ x R'^ G R^ with SE(3) action given by 

(i?,r)(x,y) = (i?x,i?y + r). 

If the measurements yield the sequence (x^^iysj {i = 1,...,/) this suggests that up to 
measurement error the configuration g{s) G SE{3) satisfies 

^(Si)(e3(0),0) = (x,^,y,J. 

The task of modelling the configuration g{s) can then be cast in the form of an inexact 
trajectory planning problem in the left-action, left-reduction form with cost functional 

%,^°,^\ A^°,^^] = £' ^ 11^^ - adt„(^° - z)||% (^°,r^^ -O + (^Sf - ^ V' 

1 ' 
+ ^Ell^(^*)(^3(0),0)-(x,,,y,jf . 

i=l 

Remark 3.3. Due to the anisotropy in velocity space, expressed by the vector z, the Euler- 
Poincare equation (3.12) is not the reduced geodesic equation for the curve g{t) with respect 
to the metric defined by K . Consequently, the solution to the inexact trajectory matching 
problem is not a Riemannian cubic spline. 

4 Geometric discretization 

The purpose of this section is to illustrate how the higher-order Hamilton-Pontryagin prin- 
ciple offers a direct route towards geometric numerical integrators. All one needs to do, in 
essence, is to provide a geometric discretization of the constraint TRg~ig — .^^ = and define 
a discrete Hamilton-Pontryagin principle accordingly. For first order variational problems 
this idea was introduced in [BRM09], building on the general theory of variational integrators 
(see [MWOl] for an extensive review). We follow in the footsteps of [BRM09] to treat second 
order problems. Third and higher orders can be dealt with in a similar fashion. 

Our main motivation for the development of a geometric integrator lies with the exact 
momentum behaviour exhibited by the discrete solution curves. Not only does the momentum 
conservation on open intervals (2.24) translate exactly into the discrete picture, but the 
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behavior at the nodes is given by discrete versions of the continuous time node equations 
(2.10)-(2.13). As a consequence, one can obtain discrete analogues of Theorem 2.2 and 
Corollary 2.3. This means that the numerical search for the optimal initial value of the 
momentum /i° can be restricted to a linear subspace of Q* of the same dimension as the data 
manifold Q. As we shall see below, the variational nature of the integrator also means that 
the discrete flow map preserves the canonical symplectic form on T*{TG). 

4.1 A geometric integrator 

In discrete mechanics the time axis [tQ,ti] is replaced by a set of discrete time points tk = 
to+kh, k = 0, . . . ,N, where h is the step size and t; = tQ+Nh. We use integers iVj, i = 1, . . . ,1, 
as node indices tj = to + Nih. For convenience we also define Nq := 0. We will need a map 
T : Q -^ G that approximates the Lie exponential and is an analytic diffeomorphism in a 
neighbourhood of with r(0) = e as well as T(^)r(— ^) = e for all 7763. An example 
is the Cayley transform, which is a second-order approximation of the Lie exponential in 
quadratic matrix Lie groups. This includes the applications discussed in Section 3. The 
Cayley transform is defined as 

r(0 = (e-e/2)-^(e + e/2). 

This and further examples are discussed in [BRM09]. 

Since r is an approximate of the Lie exponential, a simple way of discretizing the con- 
straint TRg-ig — ^° = is to require that gk+i = T{hC,^)gk, where h is the size of a time step. 
Similarly one may translate ^^ = ^^ to ,^2+i = ^2 + ^^a-- With these considerations in mind 
we define a discretized version of the action functional (2.2) on discrete path space C^. Let 

Cd = {go,eo,eo, {9k,ek,ek,f^ll^l)k=i} = (^ x 2s) x (G x 2s x 2s*)^, 
then we define 5*^ : C(i — > R as 



Sd = h 



'N-l 



k=0 

I 



E ^(^°' ^^•) + ( /^°+i' r"'(^fc+i^. ') - 4° ) + (f^Ui, k^k+i - 4°) - e^ 



+ ^2Y.d'(9N.Qo,Qu). (4.1) 

In analogy to continuous time we assume that go = e and ^q ^i's given. 

The discrete Euler-Lagrange equations follow from Hamilton's principle SS^ = 0. That is, 
they characterize paths 7 e Cd, for which SS^ := ^7(5') = for all variations ^7 G T^Cd with 
SQo = and 5^q = 0. In the process of computing SSd we need to calculate 6T^^{gk+ig^^)- 
For that purpose it is convenient to introduce the left-trivialized differential of r at ^ G s, 

whose inverse we denote by Dt^^. By taking a derivative of t{^)t{—^) = e one can show 
that [BRM09] 

Dr^^ = DtzI o Ad^(^) . (4.2) 
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Denoting r]k '■= {^gkjQk ^ '^^ ^^^^ 

= ^r{fe50)r-^(r(/i^°) Ad^(_;,^0) r]k+i) - T^(h^o)T-\r{h^'^^)r]k] 



h^i 



hei 



Dr_l^o{rik+i) - Dr^^oiVk), 



where in the last equality we used (4.2). Introducing the quantities 



/^o '■= ^} + ^/i? - ^TTo ' 

"SO 



^° := (Dr-loTfil 



and 



fil:={Dr:' r/i° {k = l,...,N), 



(4.3) 



(4.4) 



we obtain, after rearranging terms, 

'N-l 



SSd = h 



Si 



1.1 1 



E/ -0 -■- 1 -■- 1 rj^O \ / "^ 1 Tfl 

( TTo ~ /^fc+1 + r/^fc ~ "r/^fc+i' '^sfc / + \ 771 ~ l^k+iJ ^sfc 



1 



6i 



sa 



+ ( 1/^° - l-iDr,lnDr_,^or^^l^, ) + (5/i°+i, . . .) + (5/ii+„ . . .) 



+ h 



^ - f^l s^o) + /^ (5/ii, ■■■) + h (5/ii, . . .) + (/ijv, 5ej^> - if^l S^o) 



1 

i=l 



where we used abbreviations d^^ = d{gNiQo,Qu) and di^Ar^ := did{gN.Qo,QtJ. The Euler- 
Lagrange equations are therefore composed of the following equations. The constraints 



gk+i = r{Hl)gk. 4°+i = ^1 + ^1 (A; = O, . . . , iV - l), 
the discrete equations for the Legendre-Ostrogradsky momenta 



^A+l = lA- ^(^^-h««)>°+i + h 



se. 



{k = l,...,N-l] 



5i 



^-^^+1 = (fc = 0,...,iV-l), 



the discrete version of the Euler-Poincare equation for interior indices k ^ Ni {i = 1, 



and for node indices k = Ni (i = 1, ...,/ — 1) 



/.+! = iDr:l^onDT,^or (^/x° + ^J<3(di4)) . 



(4.5) 

(4.6) 

(4.7) 

..,/) 
(4.8) 

(4.9) 
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Finally, 

fi% + ^J^{d,d^) = 0, (4.10) 



a 



/ijv = 0. (4.11) 

A solution 7 G C^ to (4.5)-(4.9) is said to have initial conditions {go,^o,f^Q, fJ'o) G Gxqx2q* = 
T*{TG), using the definitions (4.3). If the Lagrangian i is hyperregular, then (4.7) can be 
solved for C,l- This means that ^l can be eliminated from equations (4.5)-(4.9), which can 
subsequently be integrated for given initial conditions. If in addition equations (4.10) and 
(4.11) are satisfied, then 7 is a critical point of the action functional 5*^. 

Remark 4.1. In Remark 2.1 we mentioned that besides the above left-action, right-reduction 
case, three other cases were available to be considered. In a similar manner to what we 
observed in that remark the right-action, right-reduction case introduces changes to equations 
(4.9) and (4.10). These become 

f^Li = {DrZl^,^*{Dr^^.y [^^l - ^ Ad;.. J«(di4)) and 
/^-^Ad* .J«(dit/^)=0, 

respectively. The remaining two cases (left-action, left-reduction; and right-action, left- 
reduction) are equivalent to the first two in just the same way as explained in Remark 2.1. 

4.2 Geometric properties 

As in Section 2.3 the interior equations are conveniently analyzed by omitting the mismatch 
penalty term (4.1) from the action functional, so that 



Jd = h 



E ^(^°'^^) + ('^°+i' {-^'^au+ig^") - 4°) + (/^Ui' ^(4%i - 4°) - e. 



The arguments that surround equation (2.23) and show symplecticity of the continuous time 
fiow can then be applied in a straightforward manner to the discrete case. Indeed, interior 
equations (4.5)-(4.8) define a fiow map F^ : G x g x 20* — )■ G x g x 20* , which integrates 
a solution 7 for given initial conditions. That is, {FdY^gQ.^Q, n^, ii^) = ((7fc,^°,/i°,/i^) or 
more succinctly {Fd)^{'~fo) = Ik- We restrict J'd to solutions of (4.5)-(4.8) and express it as 
a function J7^ initial • T*(TG) — )■ R of initial conditions 70 G T*{TG). This means that if 
7 G Cd is the solution obtained by integrating 70 then ^^^ initial (^0) ~ J'd{l)- It follows that 

d:/,,initial(^7o) = (/^jv,^^)^) " {^^lK) + {i^I.Vn) - W,m) = {{{FafTO - e){5^,\ 

where 9 is the canonical one-form on T*{TG). Taking an exterior derivative shows that the 
canonical symplectic form w = d^ is preserved by F^. Hence the discrete fiow map F^ is 
symplectic. 
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Similarly the observations given in Section 2.4 can be translated to the discrete picture. 
Indeed we pointed out in the paragraph of equation (2.25) how to obtain the conservation of 
the momentum Jk = Ad* /i° from a variational perspective. These arguments can be applied 
to the discrete variational principle to show that Ad* ^ f^t+i ~ -^^Ik f^k ^*^^ interior indices 
k 7^ Ni. Alternatively, a manipulation of equation (4.8) using (4.2) shows that 

/^^i = iDr:l^,^*iDr,,orf^l = Ad;(_,,0) /i° = Ad;^-.^ /.° = Ad;-^^ Ad; ^^l 

and therefore Ad*^^^ /i^+i = Ad*^. /i^. 

Remark 4.2. The symplectic flow equations on open intervals can equivalently he discretized 
in the purely Lagrangian picture by following [CJdD12j. One chooses a suitable discrete 
Lagrangian L^iGxCxG— >R and applies variational calculus to the discrete action sum 

N-2 

Sd = y ^ Lci{gk,gk+i,gk+2)- 

fc=0 

Let us define ^ : G x G ^ Q by ^((71,(72) = h^^r^^ {g2gi^) o,nd set 

Ld{gk,gk+i,gk+2) ■= hi{^{gk, gk+i), h'\^{gk+i, gk+2) -^{gk,gk+i)))- 

Boundary conditions being equal, the resulting optimal curve {go, . . . ,gN) is the same as the 
one we obtained from the discrete Hamilton-Pontryagin principle. The Hamilton-Pontryagin 
principle has an advantage in situations where more sophisticated discretizations of the con- 
straint TRg-ig = ^^ are chosen. For example, in the preliminary study [BHMJl] Runge- 
Kutta-Munthe-Kaas methods were used to introduce a class of such integrators. Those in- 
tegrators can still be understood in the purely Lagrangian framework, however the definition 
of the corresponding function Lii{gk,gk+i,gk+2) ^s implicit in that evaluating it requires solv- 
ing a variational problem. The Hamilton-Pontryagin approach circumvents this difficulty by 
building the discretization of the constraint into the variational principle from the outset. 

The node equation (4.9) reflects in a geometrically consistent way the jump discontinuities 
of Ad*/i° given in (2.29). Indeed, (4.9) says that 



Jfc+i = Jfc + ^Ad* J«(dirf, 



k 



when k = N^. Moreover, the final time conditions (4.10) and (4.11) are exact analogues of 
(2.12) and (2.13). See Figure 4.1 for an example. Putting everything together leads to discrete 
versions of Theorem 2.2 and Corollary 2.3. The proofs are analogous to the continuous time 
case. 



Theorem 4.3. 


Fork = 


0,. 


..,iV 




/. 


= - 


4Ad*_ 


Corollary 4.4. 


Fork = 


= 0, 


...,N 



^U<^,rf;v,Ad;^^ J«(did^^ 



for all p G Qg^Q^. 
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Fig. 4.1: Momentum norms. For the interpolating discrete cubic of Figure 3.2, the plot shows the norms 
of the momenta ^1 and /ij,. The norm of /i^ displays the momentum discontinuities at node indices as well 
as exact conservation at interior indices, in accordance with (4.8) and (4.9). The norm of ^\ demonstrates 
continuity, as found in (4.G). Both curves respect terminal conditions (4.10) and (4.11). 



4.3 Practical remarks 

The integrator derived above provides a way of finding a numerical solution to the inexact 
trajectory planning problem. 

The discrete equations of motion (4.5)-(4.9) can be employed to express the action func- 
tional J'd as a function i7^ initial • T*(TG) — ;■ R of initial conditions 70 = {go,^Q, fJ^Q, fJ^l) G 
T*(TG) = G X X 20*. The minimizer in the space of initial conditions can then be deter- 
mined by a gradient descent method. Since go = e and ,^g are given, the minimization is in 
effect only over the variables (/Xq, /Xq). By Corollary 4.4 the optimal yUg lies in the subspace of 
Q* that annihilates Qq^, to which the search can therefore be restricted. On the other hand 
one still needs to consider all of q* for the optimization of yuj. 

The gradient of i7^ initial "^^^ ^^ estimated via finite-difference methods. However, this 
requires the repeated forward integration of (4.5)-(4.9). The number of such integrations 
increases with the number of dimensions of the Lie group G, and for higher dimensional 
systems this quickly becomes unfeasible. Such difficulties can be circumvented using the 
standard method of adjoint equations, which can be easily implemented for the geometric 
discretization presented here (a detailed derivation is provided in the Appendix). Then the 
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exact gradient is obtained at the cost of integrating twice (once forward and once backward) 
a system of equations of the same complexity as the forward equations. 

The significance of the exact preservation of final time constraints (2.12) and (2.13) in the 
form of (4.10), (4.11) is to provide verification that a (local) minimum has been found. When 
the tolerance is tight (a is small) and in the absence of a good initial guess it may occur 
that the algorithm tends to a local minimum rather than the global one. Suitable initial 
guesses can be computed using a homotopy strategy. This means a step-by-step reduction 
of a, where the optimum at one value of a is taken as the initial guess at the next smaller 
value. 

We used this algorithm to generate the figures in this paper. For simulations of quantum 
splines (see Section 3.3) using the same methods we refer to [BHM12]. 



5 Discussion and Outlook 

Discussion. This paper has discussed a type of inexact trajectory planning problem whose 
optimal curves are required to pass near a sequence of fixed target positions at designated 
times within a certain tolerance. In Section 2 a new derivation of the Euler-Lagrange equa- 
tions for this type of problem was obtained by applying reduction by symmetry to a higher- 
order Hamilton-Pontryagin principle. This approach provided a new geometric interpretation 
of the previously known node equations in terms of Legendre-Ostrogradsky momenta. The 
highest order momentum was seen to undergo discontinuous jumps at the node times as a 
consequence of a partially broken Lie group symmetry. This was the content of Theorem 
2.2 and Corollary 2.3. In Section 3 several applications of the theory were discussed, which 
summoned the inexact trajectory planning problem both from a control theoretic viewpoint 
(quantum splines) as well as in the context of a type of inverse problem (rigid body splines, 
macromolecular configurations). Finally, Section 4 was concerned with the numerical ap- 
proach to solving the problem at hand. The reduced Hamilton-Pontryagin principle was 
taken as the starting point and a geometric discretization of the Euler-Lagrange equations 
was obtained, which led to exact momentum behavior and discrete versions of both Theorem 
2.2 and Corollary 2.3. This meant in particular that the search for the optimal initial value of 
the highest order momentum could be restricted to a subspace of the dual of the Lie algebra 
of the Lie group G, whose action describes the motion. 

Outlook. This work invites further development in several directions. First, one can show 
that the discrete flow map derived in Section 4 is accurate only to first order in step size 
h. The development of geometric methods with a higher degree of accuracy would be desir- 
able. The class of integrators presented in the preliminary study [BHIMll] may prove useful 
in this regard. An alternative possibility is the Lagrangian approach of [C.JdD12] together 
with suitably exact approximations of the Lagrangian function. Second, the final step in the 
numerical optimization of the cost functional was based on a shooting method with gradient 
descent in the space of initial conditions. A comparison with more sophisticated methods of 
nonlinear programming would be a useful guide for further development. Third, the example 
of the molecular strand (Section 3.4) was solved as a problem of statics. Adding the consid- 
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eration of dynamics of the strand brings one into the realm of so-called G-Strands [HIP 12], 
for which the inexact trajectory planning problem may prove interesting. Finally, for appli- 
cations in computational anatomy one must deal with infinite dimensional diffeomorphism 
groups. Extending the framework of this paper to infinite dimensions is therefore another 
challenge that lies ahead. 

Acknowledgements We thank M. Bruveris, F. Gay-Balmaz, H. O. Jacobs, M. Leok, L. 
Noakes, T. S. Ratiu, J. Vankerschaver and F.-X. Vialard for several enlightening discussions 
of this material. DDH thanks the Royal Society for a Wolf son Research Merit Award and 
the European Research Council for Advanced Grant 267382. 

A Gradient calculation via adjoint equations 

In order to implement an efficient descent method for J7^ initial ^^ ^^ useful to have an expres- 
sion for its gradient. One may obtain a gradient estimate using finite difference methods. One 
drawback lies with the inaccuracies inherent in the estimation. Moreover, if the dimension of 
the Lie algebra Q is large such estimations quickly become computationally costly. Both of 
these drawbacks can be circumvented in our case by the use of adjoint equations. In this way 
we obtain an exact expression for the gradient of J'^ initial' ^^ ^ computationally efficient way. 
We will derive the system of adjoint equations now. For simplicity we treat only the special 
case where i = | ||^^|L, where 7 denotes an inner product on q and ||.|| the corresponding 
norm. Moreover we will only treat the left-action, right-reduction case, however the others 
can be obtained in the same way. We recall from (4.5)-(4.9) the equations of motion 

g,+^ = Tihe,)9k, 4V = 4° + /^(/^Ui)* (A.l) 

^^l+l = f^l- KDr.neJ l^l+i^ (A.2) 

/x^+i = {DrZl^^YiDr^^^* (/, + A,(^,Qo)) , (A.3) 

where we introduce functions A^ : Q — t- 0* for A; = 0, . . . , A defined as 

Aiv.(g):=^J^(dirf(^iv.Qo,QO 

when k G {Ai, . . . , A;} and A^ = otherwise. 

Let us define an augmented functional Q, in which these equations are paired with La- 
grange multipliers. These Lagrange multipliers will be denoted {P^ , P^ ,V^ ,V^) G 2^* x 20 
for /c = 1, . . . , A. Let us introduce the shorthand notation x representing the discrete path 
(fl'fc) 'Cfc) /^fc) lAt)k=o ^"^^ ^ representing the ensemble of Lagrange multiplier (P°, P^, V^?, V^)^^^. 
The augmented functional Q is given by 



Q{x,X) = h 



N-l 



k=0 



+ (PfcVi, 4%i - 4° - Hk^Uif) + UUi -k^l + h{DT,,^orf,l^,, v,%. 
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1 \ 



{Dr_,^orf,l_,,-{Dr,^or{fil + A,{g,Qo)),V,\,) 



1 ' 



No constraints are assumed here, apart from the prescribed initial velocity ^q and Qq = e. It 
is clear that for any choice of Lagrange multipliers A we have Q{x,X) = iT^ initial '-^O'/^o)' ^^ 
long as X satisfies (A.1)-(A.3) for given initial values fi^, fi^. A tedious, but straightforward 
calculation shows that 



SGix, A) = -h {6fil V,') - h (6fil, Dr^^nVl 



(A.4) 

if X satisfies (A.1)-(A.3) and A is a solution of the adjoint equations. We describe these now. 
We introduce functions K^ : 20 x g* — )• g* by the defining relation 

d 

{Kl^V, p)= -^ ((Dr±;,(5+,p))>, V) , for all e, V^, P e s, /x G 0*. 

•^^ £=0 

Moreover, for A; = 0, . . . , A^ we define functions Ak : <5 x — ?• 0* by 

(AA,.(exp(£r/)g),p), 

e=0 



(A(g,p),^) 



de 



for all g G Q and p, rj E Q. The adjoint equations consist of conditions at the final time point, 

(A.5) 
(A.6) 



pO 






and the following equations for fc = 1, . . . , A^ — 1, 

P° = {Dr_^^oJ* [{Dt^IoTP^+i + AkiQkQo, Dr^^oV,\,) - h-'Ak{gkQ 

Pk = Pk+l + "'-' fc+1 ~ (,°,fJ.l ^ k+1 ~ e,l,ij!}_ , , ^k+l + -^£0,^9 

n° = ^.°M - (pD» + KPDK " 



-"fe+i 



0,/.0+A;,(3fe(50)^^+l 



-hV^ 



Dr:l,o_Dr,,A^,\,. 



(A.7) 

(A.8) 

(A.9) 
(A.IO) 



These equations are posed backwards. That is, solving the adjoint equations entails initial- 
ising the Lagrange multipliers at time point A^ according to (A.5)-(A.6) and then iterating 
backwards from k = N to k = 1 using (A.7)-(A.10). 

We now obtain an expression for the gradient of i7^ initial f^^)'^ (^•4)- Indeed, let 
(/Xo(£),/io(e)) be a variation of initial conditions (/Xo,/Xo), and let x{e) be the correspond- 
ing set of solutions to (A.1)-(A.3). Let A be a solution to the adjoint equations (A.7)-(A.10) 
for X = x(0), then 



"^d,initial 



de 



e=0 



^d,initial(/^o(^)'/^o(^)) 



de 



Qixie),\) 



e=0 



= -h {6fil, V,') - h (,5^°, Dt^^oV,') . 
From this we can read off the gradient, 

"^djnitial , r-. taI "•-'d,initial 



5/ig 



-hDT^,.oyl 



6fil 



-hV^. 
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